function h = base_plot_u(base)
% h = base_plot_u(base)
%
% Plot the primal basis of base

 meshg = load('plot_mesh.mat');
 % increase the resolution
 pgn = zeros(2,size(meshg.tg,2)*size(meshg.pg,2));
 for i = 1:size(meshg.tg,2)
   x = meshg.pg(1,meshg.tg(1:3,i));
   y = meshg.pg(2,meshg.tg(1:3,i));
   BK = [x(2)-x(1) , x(3)-x(1) ;
         y(2)-y(1) , y(3)-y(1) ];
   XY = BK*meshg.pg;
   XY(1,:) = XY(1,:) + x(1);
   XY(2,:) = XY(2,:) + y(1);
   pgn(:,(i-1)*size(meshg.pg,2)+1 : size(meshg.pg,2)*i) = XY;
 end

 U = zeros(base.pk,size(pgn,2));
 for i = 1:base.pk
   p_si = base.p_s{i};
   for j = 1:size(p_si,2) % loop over the monomials
     U(i,:) = U(i,:) + p_si(1,j)*(pgn(1,:).^p_si(2,j)) .* ...
                                 (pgn(2,:).^p_si(3,j));
   end
 end

 for i=1:base.pk
   h(i) = figure;
   hh = plot3(pgn(1,:),pgn(2,:),U(i,:),'o');
   axis equal
   title(['Basis function ',num2str(i)])
 end

return

